% This code estimates 2nd-period rate stability regulation costs by MLE
function [rs_estimate, cost2_1, cost_rs, cost2_0_mat] ...
    = rstability_estimation(s2_mkt_input, delta_p_input, ...
    in_regression_input, reg_after, fc_value, N_st, mkt_input)
%% initial values
if fc_value==1 || fc_value==2
    % fixed cost
    cost2_0_ini = ones(fc_value,1);

    % mu and sigma of log normal variable cost
    mu_cost2_ini    = 1;
    sigma_cost2_ini = 1;

    % stack
    rs_ini = [cost2_0_ini; mu_cost2_ini; sigma_cost2_ini];

elseif fc_value==N_st % baseline specification
    if isfile('rs_estimate_ini_N_st.mat')==1
        % load existing estimates
        temp = load('rs_estimate_ini_N_st.mat', 'rs_estimate');
        rs_ini = temp.rs_estimate;
    else
        rs_ini = ones(N_st+2,1);
    end
end

% bounds
LBar = -1e2*ones(size(rs_ini));
UBar = 1e2*ones(size(rs_ini));
LBar(end) = 0; 

%% estimation
% function to be minimized
ff = @(rs)rstability_likelihood(rs, s2_mkt_input, delta_p_input, ...
    in_regression_input, reg_after, fc_value, N_st, mkt_input);

options = psoptimset('Display', 'iter', ...
    'UseParallel', 'always', 'CompletePoll', 'on', 'Vectorized', 'off',...
    'TolMesh', 1e-5, ...
    'MaxFunEvals', 100000*length(rs_ini),...
    'MaxIter', 100000*length(rs_ini));

% run the optimization
disp(' ')
disp('****************')
disp('Rate stability regulation cost computation')
disp(' ')
disp('MLE running...')
tic0=tic;
[rs_estimate, fval, exitflag, output] = patternsearch(ff, rs_ini, [],[],[],[], ...
    LBar, UBar, [], options);
toc0 = toc(tic0);

%% display estimates
if fc_value==1
    vartype = {'Fixed cost',  'Log normal mu', 'Log normal sigma'};
    cost2_0     = rs_estimate(1);
    mu_cost2    = rs_estimate(2);
    sigma_cost2 = rs_estimate(3);

    cost2_0_mat = cost2_0*ones(size(s2_mkt_input)); % NN x KM

elseif fc_value==2
    vartype = {'Fixed cost, pre RSR2000', 'Fixed cost, post RSR2000', 'Log normal mu', 'Log normal sigma'};
    cost2_0_pre  = rs_estimate(1);
    cost2_0_post = rs_estimate(2);
    mu_cost2     = rs_estimate(3);
    sigma_cost2  = rs_estimate(4);

    % fixed cost as a vector, NN x 1
    cost2_0_vec = (reg_after==0)*cost2_0_pre + (reg_after==1)*cost2_0_post;
    cost2_0_mat = repmat(cost2_0_vec, 1, size(s2_mkt_input,2));

elseif fc_value==N_st % baseline specification
    cost2_0_mkt  = rs_estimate(1:end-2); % N_st x 1
    mu_cost2     = rs_estimate(end-1);
    sigma_cost2  = rs_estimate(end);

    cost2_0_vec = zeros(size(s2_mkt_input,1), 1); % NN x 1
    for ii=1:size(s2_mkt_input,1)
        mm=mkt_input(ii);
        cost2_0_vec(ii)=cost2_0_mkt(mm);
    end
    cost2_0_mat = repmat(cost2_0_vec, 1, size(s2_mkt_input,2)); % NNx KM
end

if fc_value==1 || fc_value==2
    disp(' ')
    disp('                                  Estimate    Initial       LBar       UBar ')
    for hh=1:length(vartype)
        fprintf('%-25s %15.4f %10.4f %10.0f %10.0f\n', ...
            vartype{hh}, rs_estimate(hh), rs_ini(hh), LBar(hh), UBar(hh));
    end
end

if fc_value==N_st
    % fixed cost: summary stats
    disp(' ')
    disp('Fixed cost estimates (market-specific): summary stats')
    disp(['Before RSR 2000, mean FC =', num2str(mean(cost2_0_vec(reg_after==0)))])
    disp(['After RSR 2000, mean FC =', num2str(mean(cost2_0_vec(reg_after==1)))])
    disp(' ')
    disp(['Before RSR 2000 & in event study sample, mean FC =', num2str(mean(cost2_0_vec(reg_after==0 & in_regression_input==1)))])
    disp(['After RSR 2000 & in event study sample, mean FC =', num2str(mean(cost2_0_vec(reg_after==1 & in_regression_input==1)))])

    % mu and sigma of log normal 
    disp(' ')
    disp('                                  Estimate    Initial       LBar       UBar ')
        fprintf('%-25s %15.4f %10.4f %10.0f %10.0f\n', ...
            'Log normal mu', rs_estimate(end-1), rs_ini(end-1), LBar(end-1), UBar(end-1));
        fprintf('%-25s %15.4f %10.4f %10.0f %10.0f\n', ...
            'Log normal sigma', rs_estimate(end), rs_ini(end), LBar(end), UBar(end));
end

%% compute point estimates of c1_jk
exact = s2_mkt_input./delta_p_input; % NN x K*M
exact = exact(delta_p_input>0);

threshold = (s2_mkt_input.^2)./(2*cost2_0_mat);
threshold = threshold(delta_p_input==0);

% store in NN x K*M array
cost2_1 = zeros(size(delta_p_input)); 

% when there is a price increase, then c1_jk = s_jk2/(p_jk2-p_j1)
cost2_1(delta_p_input>0) = exact;

% when there is no price increase, then c1_jk = E(c|c>(s2_mkt.^2)./(2*cost2_0_mat))
cost2_1(delta_p_input==0) = exp(mu_cost2+(sigma_cost2^2)/2)...
    *normcdf((mu_cost2+sigma_cost2^2-log(threshold))/sigma_cost2)...
    .*((1-normcdf((log(threshold)-mu_cost2)/sigma_cost2)).^(-1));

%% compute total rate stability regulation cost
cost_rs_temp = cost2_0_mat + 0.5*cost2_1.*(delta_p_input.^2);
cost_rs_temp = cost_rs_temp(delta_p_input>0);

% store in N.ist x KM array
cost_rs = zeros(size(delta_p_input)); 
cost_rs(delta_p_input>0) = cost_rs_temp;






